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ABSTRACT 

The global scales of solar convection are studied through three-dimensional simulations of compressible 
convection carried out in spherical shells of rotating fluid which extend from the base of the convection 
zone to within 15 Mm of the photosphere. Such modelling at the highest spatial resolution to date 
allows study of distinctly turbulent convection, revealing that coherent downflow structures associated 
with giant cells continue to play a significant role in maintaining the strong differential rotation 
that is achieved. These giant cells at lower latitudes exhibit prograde propagation relative to the 
mean zonal flow, or differential rotation, that they establish, and retrograde propagation of more 
isotropic structures with vortical character at mid and high latitudes. The interstices of the downflow 
networks often possess strong and compact cyclonic flows. The evolving giant-cell downflow systems 
can be partly masked by the intense smaller scales of convection driven closer to the surface, yet 
they are likely to be detectable with the helioseismic probing that is now becoming available. Indeed, 
the meandering streams and varying cellular subsurface flows revealed by helioseismology must be 
sampling contributions from the giant cells, yet it is difficult to separate out these signals from those 
attributed to the faster horizontal flows of supergranulation. To aid in such detection, we use our 
simulations to describe how the properties of giant cells may be expected to vary with depth, how 
their patterns evolve in time, and analyze the statistical features of correlations within these complex 
flow fields. 

Subject headings: convection, turbulence, Sun:interior, Sun:rotation 



1. INTRODUCTION 

The highly turbulent solar convection zone serves as 
a laboratory to guide our understanding of the complex 
transport mechanisms for heat and angular momentum 
that exist within rotating stars. One challenge is to ex- 
plain the strong differential rotation that is observed in 
the Sun, and is likely to be also realized in many other 
stars. Another concerns the Sun's evolving magnetism 
with its cyclic behavior, which must arise from dynamo 
processes operating deep within its interior. Both en- 
courage the development of theoretical models capable 
of studying the coupling of convection, magnetism, ro- 
tation and shear under nonlinear conditions. We have 
approached these challenges by turning to numerical sim- 
ulations of turbulent convection enabled by rapid ad- 
vances in supercomputing, and to helioseismology that 
provides an observational perspective of the interior dy- 
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namics. We report here on our three-dimensional simu- 
lations of compressible convection carried out in rotating 
spherical shells that capture many of the attributes of 
the solar convection zone. The evolving solutions dis- 
cussed here are obtained from the most turbulent high- 
resolution simulations conducted so far on massively par- 
allel machines. Such modelling permits us to assess, 
with hopefully increasing fidelity, the likely properties 
of large-scale convection expected to be present over a 
wide range of depths within the solar interior. In this 
paper we will describe and analyze the features of such 
giant-cell convection, and in a subsequent paper assess 
how signatures of these flows may be searched for using 
helioseismic probing. 

Helioseismology has shown that a broad variety of so- 
lar subsurface fiows are detectable in the upper reaches 
of the solar convection zone. These range from evolv- 
ing meridional circulations, to propagating bands of 
zonal flow speedup, to varying cellular flows and me- 
anderi ng streams involv ing a wide range of horizontal 



scales (iHaber et al.' 2002. 2004; Zhao fc Kosovic hev 2004: 
Hindman et"aIT ,2004. ,2006: .Kqmm et al.l 12004 2005, 
20071 iGonzalez-Hernandez et al.l 120061 ). Such detailed 
probing of flows, loosely designated as solar subsurface 
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weather (SSW), has become feasible through recent ad- 
vances in local-domain helioseismology that complement 
earlier studies of large-scale dynamics, such as infer- 
ences of the differential rotation, us ing global oscillation 
modes (e.g. [Thompson et al.ll2003D . In local helioseis- 
mology, the acoustic oscillations of the interior being 
sampled by high-resolution Doppler imaging of the so- 
lar surface can be analyzed over many localized domains 
to deduce the underlying flow fields, variously using ring- 
diagram, time-distan ce and holographic techniques (e.g. 
iGizon fc Birchl[200l . 

The horizontal resolution in such helioseismic flow 
probing, using for instance inversions of acoustic wave 
frequency splittings measured by ring analyses, can be 
of order 1° in sampling the upper few Mm just below 
the surface, and increases with depth, becoming of order 
4° at a depth of about 10 Mm. This suggests that one 
can search for explicit signatures of the largest scales 
of solar convection, or giant cells, which are a promi- 
nent feature in deep-shell simulat i ons of convection zone 
dynamics (e.g.lMiesch et al.l l2000t iBrun fc Toomrell2002l : 
iBrun et al.ll200 4h. but which are not readily evident as 
patterns in surface Doppler measurements. The presence 
of the fast and evolving flows of granulation and super- 
granulation, with combined rms horizontal flow ampli- 
tudes of order 500 m s~^, may serve to mask the an- 
ticipated weaker flows of giant cells. The helioseismic 
sampling at depths of a few Mm or greater, where the 
granular signal is likely to be sharply diminished and the 
supergranular one beginning to decrease, may thus af- 
ford unique ways to search for the largest scales of solar 
convection. 

We will here use our highest resolution, and thus most 
turbulent, spherical shell simulations of solar convection 
to assess the possible character of the giant cells. We rec- 
ognize that our solutions are at best a highly simplified 
view of the dynamics proceeding deep within the sun. 
The real sun may well possess more complex flows, or 
possibly even greater order in the form of coherent struc- 
tures, since turbulence constrained by rotation, spheric- 
ity and stratif ication can exhibit surprising behavior (e.g. 
lToomrd[2002l ). Further, our simulations here cannot yet 
deal explicitly with either the near-surface shear layer 
nor with the tachocline, concentrating instead on the 
bulk of the convection zone. However, we believe it pru- 
dent to use these models to provide some guidance and 
perspective for what may be sought with local helioseis- 
mic probing as the search for solar giant cells continues. 
The flows of SSW likely contain some signals from giant- 
cell convection over a rang e of depths (e.g. iHaber et al.l 
l2002t [Hindman et al.ll200^ . as d o power spectra of sur - 
face Doppler measurements (e.g. iHathawav et al.ll2000l ). 
We will in §3 discuss the nature of the convective struc- 
tures realized in our simulations, in §4 analyze the differ- 
ential rotation and meridional circulations that are estab- 
lished, in §5 show how the coherent downflow structures 
of the giant cells can be identifled and tracked with time, 
and in §6-7 consider the flow statistics and spectra of our 
global-scale convection. In a subsequent paper we shall 
concentrate on discussions of what may be required to 
try to resolve and possibly track the evolution of giant 
cells by helioseismic means. 

2. MODEL DESCRIPTION 



2.1. The ASH Code 

The anelastic spherical harmonic (ASH) code solves 
the three-dimensional equations of fluid motion in a 
rotating spherical shell under the anelastic approxima- 
ti on. Details on the nu merical method ca. n be found 
in IClune et all (|1999D and iBrun et al.l (I2004D an d a dis - 
cussion of the anelasti c appro x imation inlGoughl (1 19691) . 
iGlatzmai er fc GilmanI (|1981af) . [Lantz fc FanI (|1999D . and 
Micsch (2d0l). What follows is a brief summary of the 
physical model and computational algorithm. 

The anelastic equations expressing conservation of 
mass, momentum, and energy are given by 
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These equations are expressed in a spherical polar coor- 
dinate system rotating with an angular velocity of il*, 
with radius r, colatitude 9, and longitude (p. The corre- 
sponding unit vectors are f, 9, cf) and the velocity com- 
ponents are given by v = Vrf + vgO + v^<f). The density 
p, pressure P, temperature T, and specific entropy S are 
perturbations relative to a spherically-symmetric refer- 
ence state represented by p, P, T, and S. This reference 
state evolves in time, being periodically updated by the 
spherically-symmetric component of the perturbations. 
Since convective motions contribute to the force balance, 
the final term on the right-hand-side of equation ^ is 
generally nonzero. The gravitational acceleration g and 
the radiative diffusivity Kr are independent of time but 
vary with radius. 

The components of the viscous stress tensor 27 are 
given by 



V,, 



-2-pv 



and the viscous heating term is given by 
$ = 2pv 



(4) 



(5) 



In these expressions is the strain rate tensor and Sij 
is the Kronecker delta. Summation over i and j is im- 
plied in equation ([5]). The kinematic viscosity v and the 
thermal diffusivity k represent transport by unresolved, 
subgrid-scale (SGS) motions. In this paper they are as- 
sumed to be constant in space and time in order to min- 
imize diffusion in the upper convection zone, which is of 
most interest from the perspective oof helioseismology. 

The vertical vorticity ^ and the horizontal divergence 
A are defined as 

C-(Vxv).f (6) 



and 



A = V- [vee + v^(jy 



(7) 
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Equations ([T])-® are solved using a pseudospectral 
method with spherical harmonic and Chebyshev ba- 
sis functions. A second-order Adams-Bashforth/Crank- 
Nicolson technique is used to advance the solution in time 
and the mass flux is expressed in terms of poloidal and 
toroidal streamfunctions such that equation ([1]) is satis- 
fied at aU times. The ASH code is written in FORTRAN 
90 using the MPI (Message Passing Interface) library, 
and is optimized for efficient performance on scalably 
parallel computing platforms. 

2.2. Simulation Summary 

In previous papers based on ASH simulations we have 
investigated parameter sensitivities, convective struc- 
ture and transport, the maintenance of differential rota- 



namo processes ( 


Miesch et al.l l2000l lElliott et all 120001: 


Brun & Toomrel 


20021 iBrun et all 120041 iMiesch et all 


2006: Browning et al.ll20061. In this paper we focus on a 



single, representative high-resolution simulation and dis- 
cuss aspects of the flow field which can potentially be 
probed by helioseismology. 

The simulation domain extends from the base of the 
convection zone ri = 0.71i? to r2 = 0.98i?, where R is the 
solar radius. Thus, the upper boundary is about 14 Mm 
below the photosphere. Beyond r = 0.98i? the anelas- 
tic approximation begins to break down and ionization 
effects become important. The small-scale convection 
which ensues (granulation) cannot currently be resolved 
by any global model. We assume that the boundaries are 
impermeable and free of tangential stresses. 

At the lower boundary we impose a la t itudin al en- 
tropy gradient as discussed bv lMiesch et al.l ()2006l ). This 
is intended to model the thermal coupling between the 
convection zone and the radiative interior through the 
tachocline. If the tachocline is in thermal wind bal- 
ance as suggested by many theoretical and numerical 
models, then the rotational shear inferred from helio- 
seismology implies a relative latitudinal entropy varia- 
tion S/Cp ^ 4 X 10~^ where Cp is the specific heat at 
constant pressure. This corresponds to a temperature 
variation of about lOK, monotonically increasing from 
equator to pole. We implement this by setting 



Cp 



= 02^2" + ^aY^ 



(8) 



at r = ri , where YJ" is the spherical harmonic of degree 
t and order to. Here we take 02 = 2.7 x 10~^ an d 04 = 
-6.04 X 10""^. For further details see lMiesch et all (2006 ) . 
For the upper thermal boundary condition we impose a 
constant heat flux by fixing the radial entropy gradient. 

Solar values are used for the luminosity i* = 
3.846 X 10^'^ erg s~^ and the rotation rate fi* = 
2.6 X 10~^ rad s~^. The reference state, the gravita- 
tional acceleration g, and the radiative diffusion 
based on a 1-D solar structure m odel as described by 
IChristensen-Dalsgaard et al] (|1996f) . The density con- 
trast across the convection zone p(ri)/p(r2) = 132 which 
is more than three times higher than any previous simu- 
lation of global-scale solar convection. This large density 
contrast plays an important role in many aspects of the 
flow field, including the scale of the downflow network 
near the surface and the asymmetry between upflows and 



downflows. Previous simulations did not have sufficient 
spatial resolution to capture such dynamics. 
The spatial resolution Nr = 257, Ng = 1024, and 
— 2048 is higher than in any previously published 
simulation of global-scale solar convection. In our tri- 
angular truncation of the spherical harmonic series rep- 
resentation, this corresponds to a maximum degree of 
imax — 682. High resolution has enabled us to achieve 
turbulent parameter regimes that were inaccessible in 
previous simulations. We set v = 1.2x 10^^ cm^ s~^ and 
K — 4.8 X 10^^ cm^ s~^ throughout the computational do- 
main, yielding a Prandtl number Pr = v / n = 0.25. The 
velocity amplitude varies from about 250 m s~^ near the 
top of the shell to about 50 m s~^ near the bottom (Fig. 
[T5b). If we take the length scale to be the depth of the 
convection zone, D — 187 Mm, then the Reynolds num- 
ber near the top of the sheU is = UD/v ~ 400. The 
Rossby number, Ro = U/{2il^,D), varies from 0.26 near 
the top of the convection zone to 0.05 near the bottom, 
indicating a strong rotational influence on the convec- 
tive motions. However, the Rossby number based on the 
standard deviation of the vertical vorticity near the top 
of the convection zone, — 2 x 10~^, is much larger; 
Ro = ^ 4. Thus, the small-scale, intermittent 

downflows where most of the vorticity is concentrated 
are less influenced by rotation. 

3. OVERVIEW OF CONVECTIVE STRUCTURE 

Figure [1] is a representative example of the convective 
patterns achieved in the upper portion of the convec- 
tion zone. Near the top of our computational domain 
at r = 0.98i?, the structure of the convection resembles 
solar granulation but on a much larger scale; an intercon- 
nected network of strong downflow lanes surrounds a dis- 
connected distribution of broader, weaker upflows. The 
dramatic asymmetry between upflows and downflows can 
be attributed primarily to the density stratiflcation, and 
is a characteristic feature of compressible convection. As 
fluid flows upward, it diverges horizontally due to mass 
conservation. Upflows thus have a larger filling factor 
than downflows and are correspondingly less intense. 

By r = 0.95i? the downflow network begins to frag- 
ment, but isolated, intermittent downflow lanes and 
plumes remain. At low latitudes, many of the strongest 
downflow lanes have a north-south orientation. These NS 
(north-south) downflow lanes represent the dominant co- 
herent structures in the flow at low latitudes and we will 
discuss them repeatedly throughout this paper. They 
can be identified within the intricate downflow network 
near the surface but they are more prominent deeper in 
the convection zone. 

The downflow network near the surface evolves rapidly, 
with a correlation time of several days (iJ5]). Convection 
cells interact with one another and are advected, dis- 
torted, and fragmented by the rotational shear. At mid 
and high latitudes, downflows posses intense radial vor- 
ticity as demonstrated in Figure Oi, b. The sense of 
this vorticity is generally cyclonic, implying a counter- 
clockwise circulation in the northern hemisphere and a 
clockwise circulation in the southern hemisphere (more 
generallly, the vorticity vector is referred to as cyclonic 
if it has a component parallel to the rotation vector and 
anticyclonic if antiparallel). The vorticity peaks at the 
intersticies of the downflow network in localized vortex 
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Fig. 1. — Snapshots of the radial velocity Vr at (a) r = 0.98-R and (b) r = 0.95i?. Bright and dark tones indicate upfiow and downflow 
respectively as indicated by the color bar. The horizontal surfaces are displayed in a Molleweide projection which includes all 360° of 
longitude and in which lines of constant latitude are horizontal. Dashed lines indicate latitudes of 0°, ±30°, and ±60° and longitudes of 
0° and ±90° . 
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Fig. 2. — Snapshot of convective patterns near the surface. The four columns illustrate (a) radial velocity v,-, {b) radial vorticity (c) 
horizontal divergence A, and (d) temperature T. The orthographic projections in the top row correspond to r = 0.98R and the north pole 
is tilted 35° toward the observer. The three rows below show a 45° X 45° patch in latitude (10°-55°) and longitude at r = 0.98-R, 0.95iJ, 
and 0.92iJ. Color tables for each column are indicated below the 0.98il patch but scaling varies with depth. The scales indicated on the 
color bar correspond to the 0.98i? projections (upper two rows) while the scales used for the deeper layers are indicated below each image. 
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Fig. 3. — Radial energy flux as a function of radius, integrated 
over horizontal surfaces and averaged over time (62 days). Compo- 
nents include Fe (solid line), Fr (dashed), i^j. (dot-dash), and Fu 
(dot-dot-dash), all normalized by the solar flux Ft = Lt / (Anr^). 
The sum of these four components is shown as a thick solid line. 
Vertical dotted lines indicate the radial levels illustrated in Fig. |2] 



tubes which we refer to as high-latitude cyclones. Vortex 
sheets also occur in more extended downflow lanes. 

The cyclonic vorticity in downflow lanes arises from 
Coriolis forces acting on horizontally converging flows. 
Near the top of the convection zone there is a strong 
correlation between vertical velocity and horizontal di- 
vergence as demonstrated in Figure [2ji, c. This is as 
expected from mass conservation; as upflows approach 
the impenetrable boundary they diverge due to the den- 
sity stratiflcation and eventually overturn, with regions 
of horizontal convergence feeding mass into the down- 
flow lanes. Fluid parcels tend to conserve their angular 
momentum, giving rise to weak anticyclonic vorticity in 
diverging upflows and stronger cyclonic vorticity in nar- 
rower downflow lanes. Thus, the kinetic helicity of the 
flow, defined as the scalar product of the vorticity and the 
velocity, is negative throughout most of the convection 
zone, changing sign only near the base where downflows 
dive rge horizontally upon encountering th e lower bound- 
ary (jMiesch et al.ll2000l : iBrun et al.ll2004D . 

The thermal nature of the convection is evident in Fig- 
ure [2jj, d; upflows are generally warm and downflows 
cool. The more diffuse appearence of the temperture 
structure relative to the vertical velocity structure may 
be attributed to the low Prandtl number Pr = 0.25. The 
most extreme temperature variations arc cool spots as- 
sociated with the high-latitude cyclones. Global-scale 
temperature variations are also evident in Figure [2ji, in 
particular the poles are on average 6-8K warmer than the 
equator. This is associated with thermal wind balance 
of the differential rotation as discussed in []4l 

The correlation between temperature and vertical ve- 
locity gives rise to an outward enthalpy flux as illustrated 
in Figure [3l The convective enthalpy flux dominates 
the other flux components throughout most of the con- 
vection zone and peaks at r = 0.92i? where its integrated 
luminosity exceeds the solar luminosity by as much as 



70%. This is a consequence of the pronounced asymme- 
try between upflows and downflows. The greater inten- 
sity of the latter gives rise to a large downward kinetic 
energy flux Fk (cx: v'^Vr) which must be compensated for 
by the upward enthalpy flux (Fig. [3]). This has impor- 
tant consequences for 1-D solar structure models based 
on mixing length theory which generally neglect Fk and 
thus assume that the integrated convective enthalpy flux 
in the convection zone is equal to L*. 

Near the boundaries both Fe and Fk drop to zero due 
to the impenetrable boundary conditions. Flux is car- 
ried through the boundaries by radiative diffusion Fr 
and subgrid-scale (SGS) thermal diffusion F„, the latter 
of which is proportional to the radial entropy gradient 
dS /dr. Viscous heat transport is negligible and is there- 
fore omitted from Figure[3l Complete expressi ons for Fe, 
Fk, Fr, and Fu are given in IBrun et all |2004D . 

The variation of convective structure with depth 
throughout the entire convection zone is illustrated in 
FigureSl As noted with regard to Figurc[l] the downflow 
network near the surface loses its connectivity deeper 
down but isolated downflow lanes and plumes persist. 
The strongest lanes and plumes remain coherent across 
the entire convection zone, spanning approximately 190 
Mm and 4.9 density scale heights. The low-latitude NS 
(north-south) downflow lanes identified in Figure [1] are 
most prominent in the mid convection zone; near the sur- 
face they merge with the more homogeneous downflow 
network and near the base of the convection zone they 
fragment into more isolated plumes. By contrast, the 
high-latitude cyclones identified in Figure [H are largely 
confined to the upper convection zone. 

As in many turbulent flows, the enstrophy (the square 
of the vorticity vector) provides a useful means to probe 
coherent structures within the flow. Figure [H] illustrates 
the enstrophy in a square patch in the upper and mid 
convection zone. Near the surface, the high-latitude cy- 
clones dominate the enstrophy, and the vorticity is pre- 
dominantly radial (Fig. [SJi). The high spatial intermit- 
tency of these vortex structures produces some Gibbs 
ringing in the enstrophy field but this becomes neglible 
deeper in the convection zone. The enstrophy in the mid 
convection zone is dominated by vortex sheets associ- 
ated with turbulent entrainment which line the periphery 
of downfiow lanes and plumes (Fig. [5f>). Such horizon- 
tal entrainment vorticies also line the downflow network 
at r = 0.98i? but they are generally weaker than the 
vertically-aligned cyclones (Fig. [Sh) . 

4. DIFFERENTIAL ROTATION AND MERIDIONAL 
CIRCULATION 

A primary motivation behind simulations of global- 
scale convection in the solar envelope is to provide further 
insight into the maintenance of differential rotation and 
meridional circulation. These axisymmetric flow compo- 
nents play an essential role in all solar dynamo models 
and have been probed extensively by helioseismology and 
surface measurements. Although the focus of this paper 
is on the structure and evolution of global-scale convec- 
tive patterns, it is important to briefly describe the na- 
ture of the mean flows produced and maintained in our 
simulation. 

The differential rotation may be expressed in terms of 
the mean angular velocity = fJ* + {v^) /{rsmO) and 
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Fig. 4.— Radial velocity Vr at four horizontal levels (a) 0.98R, (fe) 0.92i?, (c) 0.85i?, and (d) 0.71R. The color table is as in Fig.[T] with 
the range indicated in each frame. Each image is an orthographic projection with the north pole tilted 35° toward the line of sight. The 
dotted line indicates the solar radius r = R. 



the meridional circulation may be described by a mass 
flux streamfunction ^E" defined such that 



and p (ve) 



1 a* 

r sin 6 dr 



(9) 



Angular brackets <> denote an average over longitude. 
Equation ^ applies when the divergence of the mass 
flux vanishes as required by the anelastic approximation. 
Time averages of fl and ^ are shown in Figure [6l 

The angular velocity profile is similar to the solar in- 
ternal rota tion profile inferred fr om helioseismic mea- 
surements (|Thompson et al.ll2003D . although the varia- 
tion is smaller and there is somewhat more radial shear 
within the convection zone. The mean angular velocity 
decreases by about 50 nHz (11%) from the equator to 
latitudes of 60 degrees, compared to about 90 nHz in 
the Sun. This difference may arise from viscous diffusion 
which, although lower than in previous models is still 
higher than in the Sun, or from thermal and mechanical 
coupling to the tachocline which is only crudely incorpo- 
rated in to this model throu gh our lower boundary con- 
ditions (jMiesch et al.ll200^ . For example, perhaps the 
tachocline is thinner, and the associated entropy varia- 
tion correspondingly larger, than what we have imposed 
(Sj2]). More laminar models have more viscous diffusion 
but they also have larger Reynolds stresses so many are 
able to maintain a stronger differential rotation, some 
with conical angular velocity contours as in the Sun 
(lElliott et al.|[2000l : iBnin fc Toomre|[200a [Miesch et all 
120061 ). A more complete understanding of how the highly 
turbulent solar convection zone maintains such a large 
angular velocity contrast requires further study. 

At latitudes above 30° the angular velocity increases 
by about 4-8 nHz (1-2%) just below the outer boundary 
(r = 0.95i?-0.98i?). This is reminiscent of the subsurface 
shear layer inferred from helioseismology but its sense is 
opposite; in the Sun the angular velocit y gradient is nega- 
tive f rom r = 0.95R to the photosphere (|Thompson et al.] 
l2003f) . This discrepancy likely arises from our impenetra- 
ble, stress-free, constant-flux boundary conditions at the 
outer surface of our computational domain, r — 0.98R. 
In the Sun, giant-cell convection must couple in some way 
to the supergranulation and granulation which dominates 
in the near-surface layers. Such motions cannot presently 
be resolved in a global three-dimensional simulation and 



involve physical processes such as radiative transfer and 
ionization which lie beyond the scope of our model. 

The meridional circulation is dominated by a single 
cell in each hemisphere, with poleward flow in the up- 
per convection zone and equatorward flow in the lower 
convection zone (Fig.lHl:;). At a latitude of 30°, the tran- 
sition between poleward and equatorward flows occurs 
at r 0.84-0.85 R. These cells extend from the equa- 
tor to latitudes of about 60°. The sense (poleward) and 
amplitude (15-20 m s~^), of the flow in the upper con- 
vection zone is comparable to meridional flow speeds in- 
ferred from local helioseismology arid surface measure- 
ments (iKomm et al."1993; 'Hathawav"1996'; 'Br aun fc Fal] 
1998; Habcr et al. 2002; Zhao & Kosovich_ev| |2004; 
Gonzalez- Hernandez et al.l 120061 ). The equatorward flow 
in the lower convection zone peaks at r ~ 0.75i? with an 
amplitude of 5-10 m s^^. 

Near the upper and lower boundaries there are thin 
counter cells where the latitudinal velocity (ve) reverses. 
The presence of these cells is likely sensitive to the bound- 
ary conditions and must therefore be interpreted with 
care. Global-scale convection in the Sun couples to the 
underlying radiative interior via the tacholine and to 
the overlying photospheric convection (granulation, su- 
pergranulation) in complex ways which are not yet well 
understood. The sense and amplitude of the meridional 
circulation is coupled to the differential rotation by the 
requirement that the time-averaged angular momentum 
transport by advection balance that due to Reynolds 
stresses. The weak counter cell near the upper bound- 
ary (where the flow is equatorward) is thus related to the 
positive radial angular velocity gradient at high latitudes 
seen in Figure [6l> and may arise from a misrepresenta- 
tion of the Reynolds stresses at the boundary. Likewise, 
the counter cell near the lower boundary may be sensi- 
tive to the absence of a tachocline and overshoot region. 
Previous simulations which include convective penetra- 
tion tend to exhibit equatorward meridional circulations 
thro ughout the lower co nvection zone and overshoot re- 
gion (jMiesch et al.ll2000l ) . Further work is needed to clar- 
ify the complex dynamics at the top and the bottom of 
the solar convection zone and what effect it has on mean 
flow patterns. 

Figure [7] illustrates angular momentum transport in 
our simulation including contributions from Reynolds 
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Fig. 5. — The enstrophy (oj^, where a; = VXv) shown for a 45° 
X 45° patch in latitude (10°-55°) and longitude at (a) r = 0.98R 
and (6) r = 0.85R. The color table is as in Fig. [T] but here scaled 
logarithmically. Ranges shown are (a) 10^^^ to 10~^ s~-^ and (6) 
10-1^ to 10'* s-2. 



stresses (RS), meridional circulation (MC), and viscous 
diffusi on (FD). These corresponding fluxes are deflned 
as (Elliott et al. '20001: iBrun fc Toomril2002l : iBrun et all 
[2004 iylicsch 2005) 



^RS 



^'''^ = pC({vr)r + {ve)e] 



= -pi^r'^ sin^ evil , 



where 



£ = r sin0 (Sir sin 6* + (v^)) 



(10) 

(11) 
(12) 

(13) 



is the specific angular momentum and primes indicate 
that the longitudinal mean has been removed, e.g. v'^. = 

Vr - {Vr). 



The total angular momentum flux through perpendic- 
ular surfaces is obtained by integrating the various com- 
ponents as follows: 



/;(r) = / F;{r,e)r^ sine dO 
Jo 



F^{r, e)rsm0 dr 



(14) 



(15) 



where i corresponds to RS, MC, or VD. Figure [7| shows 
time averages of these integrated fluxes. 

The prograde differential rotation at the equator is 
maintained primarily by equatorward angular momen- 
tum transport induced by Reynolds stresses (Fig. [7f)). 
This transport is dominated by the NS downflow lanes 
discussed in Sj3] which represent the principal coherent 
structures at low latitudes (persistent over relatively long 
times and large horizontal and vertical scales). Coriolis- 
induced tilts in the horizontally converging flows which 
feed these downflow lanes give rise to Reynolds stresses 
which transport angular momentum toward the equator 
(see Miesch 2005, Fig. 15). 

Reynolds stresses also transport angular momentum 
inward throughout most of the convection zone (Fig.lTli). 
This is a significant departure from previous simulations 
of global-scale solar convection which have exhibited an 
outward t ransport of angular momentum by Rey nolds 
stresses ( Brun fc Toomi?l2002l: iBrun et a\.\ \2004). In- 
ward angular momentum transport by convection is a 
common feature of many mean-field models in which 
it is typicallyparameterized by means of the so- called 
A-effect fe.g. pCit chatinov & Riidiger 1993: Canuto "et al.l 
11994'; Kitc hatinov fc Riidiger 2005; R iidiger et al.i i2005'). 
However, in some models this inward transport arises 
from a velocity anisotropy such th at the standard dev i- 
ation of Vr exceeds that of fe.g. iRiidiger et aLll2005f) . 
Such is not the case in our simulation where the three ve- 
locity components are comparable in amplitude through 
most of the convection zone (see Fig. [rSb). The rever- 
sal in F^^ and F^^^ near the boundaries is associated 
with the counter cells in the meridional circulation seen 
in Figure [nt- 

Advection of angular momentum by the meridional 
circulation gives rise to poleward and outward trans- 
port, nearly balancing the Reynolds stresses, while the 
transport due to viscous diffusion is relatively small. 
This approximate balance between Reynolds stresses 
and meridional circulation with regard to angular mo- 
mentum transport is that which is expected to exist 
in the convection zone of the Sun and in other stars 
where Lorcntz forces and visco us diffusio n are negligi- 
ble d^ ssoul 1978; Zahn 1992; Ell iott et al . 2000; Rcmpe] 
l2005n Micsch 2005). The curves shown in Figure [7] do not 
sum precisely to zero, indicating that there is some evo- 
lution of the rotation profile over timescales which are 
longer than the 112-day averaging interval. 

The delicate balance between F^"^ and F*^*-^ plays an 
essential role in determining what meridional circula- 
tion patterns are achieved. In previous global simula- 
tions, viscous angular momentum transport was signif- 
icant and this balance was disrupted. Circulation pat- 
terns were generally multi-celled in latitude and radius 
(jMiesch et al.ll2000l: lEUiott et al.ll2000l : iBrun fc Toomi 
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Fig. 6. — Differential rotation, meridional circulation, and mean temperature perturbation averaged over longitude and time (58 days). 
The angular velocity shown (a) as a 2-D image and (b) as a function of radius for selected latitudes as indicated. Contour levels in (a) 
are every 10 nHz and the rotation rate of the coordinate system (414 nHz) is indicated on the color bar (black line). Contours of the 
streamfunction ^ in (c) represent streamlines of the mass flux with red (black contours) and blue (grey contours) representing clockwise 
and counter-clockwise circulation respectively. The color table saturates at >!' = ±1.08 X 10^^ g s~^. Characteristic amplitudes for (vg) 
are 20 m (poleward) at r = 0.95-R and 5 m s~^ (equatorward) at r = 0.75-R. Contour levels for the temperature perturbation (d) are 
every IK. 



run et al]l2004[ ). By contrast, the circulation pat- 
terns shown in Figure [Sf) are dominated by a single cell 
in each hemisphere. 

Not only is the viscosity lower than in most previous 
simulations, but this is the first global simulation to ex- 
tend from the base of the convection zone to r = 0.98i?, 
spanning a factor of more than 130 in density ( §2.2p . 
Furthermore, we have incorporated some aspects of the 
coupling between the tachocline and the convective en- 
velope into our model by applying a weak latitudinal 
entropy vari ation at the bottom boundary ( §2.2p . As 
discussed by iMiesch et al.l (|2006) , this entropy variation 
is transmitted throughout the domain by the convective 
heat flux and the resulting baroclinicity promotes coni- 
cal angular velocity profiles which satisfy thermal wind 
balance. Such profiles minimize diffusive angular mo- 
mentum transport in radius. 

The temperature variations associated with thermal 
wind balance are evident in Figure [5Ji. The poles are 
about 6-8K warmer than the equator on average. The 
background temperature varies from 2.2x10^ K at the 
base of the convection zone to 8.4x10'* K at the outer 
boundary so the relative latitudinal variations are small, 
3 X 10-6 to IQ-'^. 

The differential rotation profile is steady in time; in- 
stantaneous snapshots appear similar to Figure [SJi but 
with more small-scale structure and somewhat more 
asymmetry about the equator. For illustration, the am- 
plitude of the temporal variations sampled at a latitude 
of 30° relative to a two-month mean is ±10 nHz to- 
ward the top of the convection zone (~ 2%), decreasing 
to ± 5 nHz toward the base. Angular velocity varia- 
tions are larger at high latitudes where the moment arm 



(rsin^) approaches zero. The amplitude and nature of 
these variations is comparable to sola r rotational vari- 
ation s inferred from helioseismology (Tho mpson et all 
l2003l ). However, in addition to more random fluctua- 
tions, the solar rotation exhibits periodic torsional oscil- 
lations which are not realized in our simulation. These 
are associated with magnetic activity which lies beyond 
the scope of our current model (e.g. ICovas et al.l l2000l : 
Spruit 2003; Rempel 2005, 2007). 

By contrast, fluctuations in the meridional circulation 
are large relative to the temporal mean, changing sub- 
stantially over the course of one rotation period as illus- 
trated in Figure [SI Variations in the axisymmetric lati- 
tudinal velocity (vg) at a latitude of 30° reach ±60 m s~^ 
at the top of the convection zone and ±10 m s~* near the 
base, as much as 300% of the two-month mean. Instanta- 
neous circulation patterns are in general multi-celled in 
latitude and radius and asymmetric about the equator. 
Some asymmetry persists even in two-month averages 
(Fig. [6]:). Large relative variations in the meridional cir- 
culation are expected because it is weak relative to the 
other flow components so it is easily altered by fluctuat- 
ing Reynolds stresses and Coriolis forces. The volume- 
integrated kinetic energy contained in the meridional cir- 
culation is approximately an order of magnitude smaller 
than in the differential rotation and approximately two 
orders of magnitude smaller than in the convection. 

Determinations of the solar meridional circulation 
from surface measurements and helioseismic inversions 
are g enerally averaged over at least on e rotation pe- 



riod (iKomm et al.lll993l : lHath awavlll99d iBraun fc FanI 

h^, ■^^=^ i 



19981: iHaber et al.ll2(M 2004: Z hao fc Kosovichevl[200l 
Gonzalez-Hernandez et al.l i2006i ) . The time variations 
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Fig. 7. — (a) Radial and (6) latitu dina l tra nsport of angular 
momentum as expressed in equations I I10I I- II151 . averaged over a 
time interval of 112 days (four rotation periods). Shown are con- 
tributions due to Reynolds stresses (RS, solid lines), meridional 
circulation (MC, dashed lines), and viscous diffusion (VD, dotted 
lines). Zero is indicated by the dash-dotted line. 



are therefore less than in the snapshots iUustrated in Fig- 
ure [Hl/'-J but consistent with comparable running time 
averages in the simulation. However, as with the angu- 
lar velocity, systematic variations in the solar meridional 
circulation associated with the magnetic activity cycle 
are not captured in this non-magnetic simulation. 

5. IDENTIFICATION AND EVOLUTION OF COHERENT 
STRUCTURES 

In ^we described recurring convective features includ- 
ing NS downflow lanes found at low latitudes and inter- 
mittent, high-latitude cyclones. In this section we discuss 
these coherent structures in more detail and address the 
lifetime, propagation, and evolution of convective pat- 
terns. 

Figure [9] illustrates the evolution of the low-latitude 
downflow network at r = 0.98ii!. Substantial changes are 
evident even over the two-day time interval between ad- 
jacent bands. Individual convection cells typically lose 
their identity after only a few days and none are clearly 
recognizable after one rotation period. This has impor- 
tant implications for subsurface weather diagrams in- 
ferred from local helioseismology (2]); temporal sampling 
of a day or less may be necessary to reliably follow the 



evolution of flow fields associated with giant convection 
cells. 

Embedded within the more rapidly evolving downflow 
network are features that persist for a month or more. 
These are the NS downflow lanes discussed in [JSj appear- 
ing in Figure [9] as dark vertical stripes although it takes 
some scrutiny to see them amid the complex smaller 
scales. These propagate in an eastward (prograde) di- 
rection relative to the rotating coordinate system as il- 
lustrated by the white arrow. 

The presence of NS downflow lanes is more readily ap- 
parent when the divergence of the zonal velocity dv^ / d4> 
is plotted as in Figure [TOl Whereas the horizontal diver- 
gence A corresponds closely to the radial velocity pat- 
terns shown in Figure [9l the zonal component alone pref- 
erentially selects structures with a north-south orienta- 
tion. Thus, the NS downflow lanes are more prominent 
and their prograde propagation and persistence over time 
scales of at least a month are evident. The propagation 
rate varies as individual lanes continually catch up to 
others and subsequently merge. 

Using Figure [10] as a reference facilitates the detec- 
tion of NS downflow lanes within the intricate down- 
flow network of Figure [9l In other words, it is eas- 
ier to distinguish the NS downflow lanes if one knows 
where to look. Furthermore, a close comparison of Fig- 
ures [9] and \W\ reveals that the horizontal scale of the 
convective cells is somewhat smaller in the vicinity of 
the NS lanes (see, for example, the downflow lane traced 
by the arrow). This is consistent with the more gen- 
eral characteristic of turbulent compressible convection 
that downflow lanes tend to be more turbulent an d vor- 
tical t han the broader, weaker upflows jBr ummcU et al.l 
11996; Brandenburg et al. 1996: Stein fc Nordl und 1998|; 
iPorter fc Woodward. 2000; .Miesch, 20051 see also Fig. O. 
Advection of smaller-scale convection cells and vortices 
into extended NS downflow lanes is apparent in anima- 
tions of the radial velocity field. 

The longitudinal position of lanes of zonal velocity con- 
vergence in the surface layers corresponds closely to the 
position of NS downflow lanes in the mid convection zone 
where they are more prominent in the radial velocity field 
(Fig.|4]). Thus, searching for lanes of zonal convergence in 
solar subsurface weather (SSW) maps inferred from local 
helioseismology might be a promising way to detect con- 
vective structures which extend deep into the convection 
zone. However, care must be taken when interpreting 
such results. Even if were isotropic in latitude and 
longitude, the one-dimensional derivative dv^/d4> would 
still exhibit a preferred north-south orientation. Thus, 
anisotropy in dv^/d(j) should not be used naively as a 
criterion for establishing the existence of NS downflow 
lanes, but it may be used to track the propagation and 
evolution of such coherent structures if they are indeed 
present. 

The NS downflow lanes are confined to latitudes less 
than about 30°. At higher latitudes the downflow net- 
work is more isotropic in latitude and longitude and pos- 
sesses intense cyclonic vorticity (5j3]). An illustrative ex- 
ample of the evolution of mid-latitude convective pat- 
terns is shown in Figure [TTl 

As demonstrated in Figures[2]and[5l intense, vertically- 
oriented, cyclonic vortices are prevalent throughout the 
downflow network in the upper convection zone. Cen- 
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Fig. 8. — Instantaneous snapshots of the meridional circulation streamlines at five times separated by an interval of 28 days. Color tables 



and contours are as in Figure |6j: but the normalization for ^ is three times larger, ±3 X 10^ 
s~^ near the top of the convection zone. 



g s 



Peak amplitudes for vg reach 60 m 



trifugal forces can evacuate the cores of the most in- 
tense vorticies, leading to a reversal in the buoyancy 
driving which siphons fluid up from below and creates 
a new upflow within the intersticies of the downflow net- 
work. This phenomenon has been referred to as dynam- 
ical buoyancy and is a chara cteristic feature of rotat- 
ing, compressible convection (iBrandenburg et al.1 Il996t 
iBrummell etallll996t iMiesch et al.ll2000D . Theresult is 
a helical vortex tube with upflow at its center and down- 
flow around its periphery. In the vertical velocity (or 
the horizontal divergence) field these structures appear 
as small convection cells, with a horizontal extent com- 
parable to that of supergranulation, about 10-30 Mm. 
Several examples of these are indicated in Figure [TT] (A). 

The formation of one of these helical convection cells 
via dynamical buoyancy is also indicated in Figure [TT] 
(B). At a (relative) time of I'', a counter-clockwise swirl 
can be seen near one of the interstices of the downflow 
network, reflecting the presence of a vertically-oriented 
vortex tube. Such cyclonic swirl is evident throughout 
the downflow network in animations of the flow field. 
One day later, a strong downflow plume develops and 
Coriolis forces continue to amplify the cyclonic vorticity. 
By the next day, centrifugal forces have evacuated the 
vortex core and reversed the axial flow. After formation, 
such upflows may spread horizontally due to the density 
stratification or they may dissipate through interactions 
with surrounding fiows. 

The horizontal spreading of a new upflow is limited 
by interactions with adjacent convection cells and by the 
need to transport heat outward and u ltimately th rough 
the boundary, as discussed by iRast] ( 19951 l2003i ). As 
can be seen in Figure [TT] (see also Fig. Wp), the strongest 
upflows occur adjacent to the downflow lanes. As a con- 
vection cell expands horizontally, the upward flow at the 
center of the cell drops, leading to a reduction in the out- 
ward enthalpy flux. Cooling of the fluid due to the up- 
per boundary condition eventually reverses the buoyancy 
driving, thus forming a new downflow lane which bisects 
and thereby fragments the existing convection cell. This 



occurs continually in our simulation as demonstrated in 
Figure [TT] (C). Similar dynamics also occur at lower lat- 
itudes, as can be seen by careful scrutiny of Figure [1] 
Such fragmentation induced by cooling near the upper 
boundary is the principle factor in determining the size 
and lifetime of the cells which make up the downflow 
network. 

Convection cells may also be squeezed out of existence, 
or collapse, via the horizontal spreading of adjacent cells 
as illustrated in Figure [TT] (D). Shearing of convection 
cells by differential rotation also limits their lifetime and 
horizontal scale. Such processes typically occur over the 
course of several days but some convection cells can per- 
sist with little distortion for nearly a week (Fig. [TT] E). 

A quantitative measure of the lifetime and propagation 
rate of convective patterns can be obtained by consider- 
ing the autocorrelation function, acf, defined as follows: 

acf(r, i, f2t,r) = 

Sel lo^ (''j ^T4'Tt)vr{r,9,(j) — fltT,t -\- r) sin 9d9d(j) 

le^ C^rir,eA,t) sineded<l> 

(16) 

where fit is the tracking rate (expressed as an angular 
velocity) , t is the temporal lag and 6i and 6*2 specify the 
desired latitudinal band (averaged over the northern and 
southern hemispheres). Results are illustrated in Figure 

The acf is unity at r = and drops monotonically 
with increasing lag. For each value of fit one may define 
a correlation time as the time beyond which the acf drops 
below a fiducial threshold, here taken to be 0.05. We then 
define an optimal tracking rate Oc as that value of f2t 
which maximizes the correlation time, Tc- Optimal rates 
and correlation times for various latitude bands are listed 
in Table [T] Also listed for comparison is the variation of 
the mean rotation rate J7 across each latitude band for 
the radius and time interval used to compute the acf. 

The acf in Figure [T2tj, corresponds to low-latitude con- 
vective patterns near the surface, such as in Figure [5] 
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Fig. 9. — Temporal evolution of convective patterns near the equator. Each image shows the radial velocity field at r = 0.98-R for all 360° 
of longitude and a latitude range of 0°— 25° . The color table is the same as in Fig. [TJi. Snapshots at different times are stacked vertically, 
with time increasing upward. The interval between snapshots is two days. To facilitate comparison, the uppermost image illustrates the 
convection structure one rotation period (27.4 days at this latitude and radius) prior to the At =28-day image immediately below it. The 
tracking rate is the rotation rate of the coordinate system, 414 nHz, but the white arrow shows a propagation rate of 450 nHz for reference. 
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Fig. 10. — Time evolution of the zonal velocity derivative dv^/d(f> at r = 0.98H. The layout and time series corresponds directly with 
Fig.|9] with bands spanning 0°-25° in latitude and time increasing upward. The color table saturates at ±200 m s~^; bright tones denote 
convergence (dv^/dtj) > 0) and dark tones denote divergence. The white arrow corresponds to an angular velocity of 450 nHz as in Fig. [9] 
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Table 1. Correlation times and Optimal Tracking Rates 





radius 


0°-20° 


20°-40° 


40°-60° 


60°-90° 


Q (nHz) 


0.98/? 


414-421 


421-408 


408-375 


375-358 


Cc (nHz) 


0.98R 


450 


430 


380 


330 


Tc (days) 


0.98iJ 


2.6 


2.0 


2.0 


2.2 


n (nHz) 


0.85R 


429-415 


415-400 


400-372 


372-358 


Qc (nHz) 


0.85i? 


440 


424 


400 


390 


Tc (days) 


0.85R 


9.0 


8.2 


8.0 


8.0 



Here the flow is dominated by the intricate, continually 
evolving downflow network and correlation times are only 
a few days. The maximum correlation time of 2.6 days 
is achieved with a tracking rate of 450 nHz which cor- 
responds to the propagation rate of the NS downflow 
lanes as indicated by the arrow in Figure [9l These NS 
downflow lanes are the longest-lived structures within the 
downflow network and they propagate faster than the 
mean rotation rate of 414-421 nHz (Table[T]). Thus, they 
are propagating convective modes as opposed to passive 
features being advected by the differential rotation. 

The NS downflow lanes are more prominent deeper in 
the convection zone and this is reflected in the acf of 
Figure fT2b. The acf is more strongly peaked at the opti- 
mal tracking rate and the associated correlation time is 
longer, Tc = 9.0 days. At 440 nHz, ^lt is somewhat less 
at r = 0.85i? than at r = 0.98i? but it is still faster than 
the local rotation rate (Table [1]). 

The prograde propagation of NS downflow lanes can 
be attributed to the approximate local conservation 
of potential vorticity, and in this sens e they may 
be regarded as thermal Ro ssby waves (jBussd 119701 : 
iGlatzmaier &: GilmanI Il981b[ ) . NS downflow lanes are 
related to banana cells and columnar convective modes 
that occur in more laminar, more rapidly rotating, and 
more weakly stratified systems and that have been well 
studied both analyticall y and numer ically (reviewed by 
IZhang fc Schubert! 120001 : iBussd 12002^ ) . In general, their 
propagation rate depends on the rotation rate, the strat- 
ification, and the geometry of the shell. 

At mid latitudes, the correlation times are somewhat 
smaller and the optimal tracking rates are slower, com- 
parable to the local rotation rate (Table [J) . Near the 
poles il and fit become less reliable because the small 
momentum arm induces large temporal variations in an- 
gular velocity. Linear theory indicates that polar convec- 
tive modes should propagat e slowly retrograde () GilmanI 
119751: iBusse fc Cuon3ll977D but it is uncertain whether 
such linear modes persist in this highly nonlinear pa- 
rameter regime. Correlation times at high latitudes are 
comparable to those at mid latitudes, about two days 
for the downfiow network in the upper convection zone 
and about 8 days for the larger-scale flows in the mid 
convection zone. 

We emphasize that statistical measures such as Tc can 
drastically underestimate the lifetime of coherent struc- 
tures within a turbulent flow such as this. It is evident 
from Figures [9l and fTOl that some NS downflow lanes per- 
sist for weeks and even months. Likewise, some higher- 
latitude convective cells at r = 0.98i? can persist for up 



to a week (e.g. Fig. [TTl E). 

6. HORIZONTAL SPECTRA AND LENGTH SCALES 

In S}3]we discussed the convective patterns realized in 
our simulations and in [JS] we described their temporal 
evolution. In this and the following section we consider 
univariate and bivariate statistics in order to gain further 
insight into the nature of the convective flows. Through- 
out this analysis we will be concerned solely with fluc- 
tuating quantities (m > 0), indicated by primes. The 
structure and evolution of mean flows is discussed in 52) 
Furthermore, we focus on the upper portion of the con- 
vection zone which is most relevent to local helioseismol- 

ogy- 

Figure [T3] shows the spherical harmonic spectra of the 
velocity and temperature fields at two levels in the up- 
per convection zone as a function of spherical harmonic 
degree £, which may be regarded as the total horizontal 
wavenumber. At r = 0.98i?, the radial velocity spectrum 
(Fig.[T5b) increases with £ approximately as reaching 
a maximum at £ ^ 80. Afterward it drops with a best-fit 
exponent of n « —5, where the power P{£) (x 

The spherical harmonic degree ^ = 80 corresponds to a 
horizontal scale L of 54 Mm, which may be regarded as 
a characteristic scale of the downflow network illustrated 
in Figures [T]i and |4ji and in the upper row of Figure 
[^h- However, a single characteristic scale is somewhat 
misleading as the downflow network exhibits structure on 
a vast range of scales. A look at the convective patterns 
in FigurefTTl for example, reveals convection cells 10°-20° 
(100-200 Mm) across as well as cyclonic vortices spanning 
only a few degrees (10-30 Mm). Meanwhile, NS downflow 
lanes can extend to latitudes of ±25° or more, spanning 
more than 500 Mm (e.g. Fig. fTO)) . 

Deeper in the convection zone, the convective scales 
are generally larger. The radial velocity spectrum at r = 
0.92 peaks a,t £ — 26, corresponding to a horizontal scale 
of about 150 Mm (Fig. [T3b. dotted line). Furthermore, 
the high-f dropoff in power is somewhat steeper than at 
r = 0.98i?, with an exponential providing a better fit 
than a polynomial; P{£) cx exp {a£), with a = —0.015. 

The horizontal velocity spectra peak a,t £ = 10-20 {L ^ 
200-400 Mm) for both radial levels sampled in FigurefTSb. 
c. These scales are larger than for the radial velocity as a 
consequence of mass conservations. Near the impenetra- 
ble boundary, the radial velocity is highly correlated with 
the horizontal divergence such that ~ A ~ £vh where 
Vh is the horizontal velocity. The Vg and v'^ spectra are 
thus somewhat flatter than v'^ at high £ {n ^ —4.6). At 
r — 0.92R, the v'g and v'^ spectra, like the v'^ spectrum, 
are best fit by exponentials with a ~ —0.016. 

The correspondence between the radial velocity Vr and 
the horizontal divergence A near the upper boundary is 
apparent when comparing frames a and d of Figure 1131 
At r = 0.98i? the two spectra are nearly identical when 
normalized by their maximum value. By r — 0.92 differ- 
ences become significant, with the A' spectrum shifted 
toward higher wavenumber relative to the v'^ spectrum. 

The vertical vorticity spectra in Figure 113b peak at 
even higher wavenumber, £ — 140 (L ~ 30 Mm) at 
r = 0.98. This reflects the presence of the high-latitude 
cyclones discussed in (!j3]) which are highly intermittent 
in space and time. Beyond this maximum, the spec- 
trum decays approximately as £~'^. Deeper in the con- 
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Fig. 11. — Temporal evolution of mid-latitude convection patterns. The radial velocity field is shown at r = 0.98iJ over a square 30° X 30° 
patch in latitude and longitude, spanning latitudes of 30°-60°. The color table is as in Fig. [T}i. The series shown spans one week, with 
an interval of one day between successive images. The patch is tracked at an angular velocity of 410 nHz. Labels indicate (A) cyclonic 
vorticies, (B) dynamical buoyancy, (C) fragmentation and (D) collapse of convection cells, and (E) the persistence of some cells for multiple 
days. 
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Fig. 12. — The acf is shown as a function of tracking rate fit 
and lag r for latitude band 0°-20° at (a) r = 0.98R and (b) r = 
0.85-??. For each tracking rate, the acf drops from a value of unity 
(yellow) to zero (white) with contour levels spaced at intervals of 
0.05. Dotted lines indicate the the optimal tracking rate and the 
associated correlation time. 



vection zone at r — 0.92R, the peak shifts toward lower 
wavenumber 100, L 40 Mm) and the spectrum 

steepens (n ^ —4). 

The spectrum of temperature fluctuations is flatter 
than the velocity field at low £, with a broad maximum 
at £ ~ 16 (260 Mm). This reflects the low Prandtl num- 
ber and the large-scale thermal variations associated with 
thermal wind balance (Q. The slope at high £ is com- 
parable to the velocity field, with n ~ —4.2 at r = 0.98 
and a 0.019 at r = 0.92i?. 

7. PROBABILITY DENSITY FUNCTIONS AND MOMENTS 

More detailed information about the nature of the flow 
field may be obtained from probability density functions 
(pdfs) as shown in Figure [TH These are normalized his- 
tograms on horizontal surfaces which take into account 
the spherical geometry. The shape of a pdf f{x) may be 
described through its moments of order n, defined as 

= J {x- {x)f f{x)dx 

= — / / ix - {x)f sixiO de d4> . (17) 
47r7o Jo 



We then define the standard deviation cr, the skewness 
5, and the kurtosis K. as follows: a = (A^^)^/^, S — 
and K. = Resuhs are illustrated in 

Figure [T51 as a function of radius. 

Near the top of the convection zone, the radial veloc- 
ity pdf exhibits a bimodal structure, with two distinct 
maxima at positive and negative Vr (Fig. 114b . solid line). 
These maxima suggest characteristic velocity scales of 
^ 6 m s~^ for upflows and ^ 30 m s^^ for downflows. 
However, these values are substantially smaller than the 
standard deviation (rms value) of v'^ which peaks at 160 
m at r = 0.96i? and then drops to zero at the upper 
boundary (Fig. [Tsh). 

The larger amplitude of the positive peak reflects the 
larger flUing factor of upflows relative to downflows. By 
r — 0.92R, the negative peak has largely disappeared 
and the negative tail of the pdf becomes nearly exponen- 
tial (dotted line). This signifies turbulent entrainment, 
whereby much of the momentum of downfiow lanes and 
plumes is transferred to the surrounding fiuid and dis- 
persed. The asymmetry between narrow, stronger down- 
fiows and broader, weaker upfiows is a consquence of the 
density stratification (Q and is manifested as a large 
negative skewness which persists throughout the convec- 
tion zone (Fig. [TSb). 

The kurtosis /C is generally regarded as a measure of 
spatial intermittency but large values can also arise from 
bimodality. A unimodal Gaussian distribution yields 
!C — 3 whereas an exponential distribution yields /C = 6. 
The v'^ pdf has an even larger kurtosis ranging from 3- 
12 across the convection zone (Fig. ITSb). reflecting both 
intermittency and bimodality. 

By comparison, the horizontal velocity pdfs shown in 
Figure fT4b . c appear more symmetric with nearly expo- 
nential tails (/C — 3-5; Fig. [TSb). The positive skewness 
of the v'^ pdf is a signature of the NS downflow lanes dis- 
cussed in !j3] and SJH Their north-south orientation and 
prograde propagation implies converging zonal flows in 
which the eastward velocities on the trailing edge of the 
lanes are somewhat faster on average than the westward 
velocities on the leading edge. 

Velocity amplitudes increase with radius due to the 
density stratification, with horizontal velocity scales 
reaching over 200 m s^^ near the surface (Fig. [T5b). 
In the mid convection zone the velocity field is nearly 
isotropic, with a characteristic amplitude of about 100 
m s^^ for all three components. 

The horizontal divergence pdf shown in Figure I14| j 
is nearly identical to the v'^. pdf in Figure [Hb at 
r = 0.98, but as with the spectra in Figure [13b . this 
correspondence breaks down by r = 0.92i?. In the 
mid-convection zone the A pdf becomes more sym- 
metric with nearly exponential tails {S — -0.13, /C — 
9.2 at r = 0.92i?). Non-Gaussian behavior such that 
/C > 3 for velocity differences and derivatives is a 
well-known feat ure found in a wide variety of turbu- 
lent flows (e.g. 'Chen et al.' '1989'; 'Castaing et al.' '19901; 
She 1991; Kailasnath ct al. 19 92; Mies ch et al. 1991; 
Jung'fc Swinnevll2005HBruno fc C"arbonelf2005h . In par- 
ticular, the pdf of velocity differences between two points 
separated in space is often modeled using stretched ex- 
ponentials f{x) oc exp(— jxp), where (3 approaches unity 
for small spatial separations (as sampled by derivatives) 
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Fig. 13. — Power spectra plotted as a function of spherical harmonic degree £, summed over all orders m > and averaged over a time 
interval of 62 days. Solid and dotted lines sample spherical surfaces r = 0.98i? and r = 0.92i? respectively. Quantities include (a-c) the 
velocity components, (d) the horizontal divergence, (e) the vertical vorticity, and (/) the temperature fluctuations. 



and becomes more Gaussian {(3 « 2) as the spatial sepa- 
ration increases. 

The vertical vorticity pdfs shown in Figure 114b also 
exhibit nearly exponential tails. However, near the sur- 
face (r — 0.98i?) the distribution is bimodal with promi- 
nent tails signifying an abundance of extreme events 
(/C = 180). These tails arise from the intense, inter- 
mittent cyclones which develop at the intersticies of the 
downflow network at mid and high latitudes as discussed 
in iJ31 By r = 0.92R, the bimodality is absent, although 
the pdf is still highly intermittent (/C = 25). This is 
consistent with Figure ^ which suggests that the high- 
latitude cyclones are confined to the outer few percent of 
the convection zone (r > 0.95i?). 

The signature of high-latitude cyclones is also present 
in the pdf of temperature fluctuations as a prominent 
exponential tail on the negative side at r = 0.98i? (Fig. 
[T4)f ). As with the radial velocity pdf, the bimodality 
disappears by r = 0.92i? due to entrainment and the 
negative tail becomes unimodal and exponential. The 
asymmetric shape of the tempature pdfs arises from the 
asymmetric nature of the convection noted in ij3l cool 
downflows are generally less space-filling and more in- 
tense than warm upflows. The temperature pdfs remain 
asymmetric {S < 0) and intermittent (/C > 3) through- 
out the convection zone, becoming most extreme near 
the surface where S — —1.6 and K. = 12. The standard 



deviation of the temperature fluctuations ranges from 0.4 
K in the lower convection zone to a maximum of 4 K at 
r = OMR. 

Correlations between vertical velocity and temperature 
fluctuations may be investigated further by means of two- 
dimensional pdfs (histograms) as illustrated in Figure 
[TBI for r — 0.98i?. Although warmer and cooler tem- 
peratures are associated with upflows and downflows re- 
spectively, the relationship is not linear. Upflows exhibit 
a prominent maximum at v'^. 20 m s^^ and T' ~ 2 
K, whereas downflows are more distributed, both in the 
range of velocity amplitudes and in the spread of tem- 
perature variations for a given u^. This spread increases 
somewhat toward higher latitudes due to the prepon- 
derence of intermittent cyclonic plumes but the aver- 
age correlation shown in Figure 1161 ! is insensitive to lati- 
tude. The reversal in the sense of the temperature vari- 
ation at high radial velocity amplitudes is due in part 
to poor statistics (few events) but it does have physical 
implications. As noted in fJSj the strongest upflows oc- 
cur adjacent to downflow lanes. Cool regions associated 
with downflow lanes tend to be more diffuse than the 
lanes themselves as a result of the low Prandtl number 
{Pr ~ 0.25). Thus, the fastest upflows can be relatively 
cool. Similarly, the fastest downflows occur in localized 
regions adjacent to warmer upflows such that the temper- 
ature fluctuations are diminished by thermal diffusion. 
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Fig. 14. — Probability density functions (pdfs) are shown for (a) v'^, (b) Vg, (c) v'^, (d) A', (e) and (/) T'. The pdfs shown correspond 
to spherical surfaces at r = 0.98R (solid lines) and r = 0.92 (dotted lines) and are averaged over time (62 days). Moments of the fluctuating 
velocity pdfs as a function of depth are shown in Fig. 1151 




Fig. 15. — The (a) standard deviation, (b) skewness, and (c) kurtosis of the three velocity components are shown as a function of radius, 
averaged over a time interval of 62 days. The solid, dotted, and dashed lines represent dJ., v'g and v'^ respectively. Horizontal dashed lines 
show skewness and kurtosis values for a Gaussian distribution (5 = 0, /C = 3) for comparison. 



Correlations between the horizontal velocity compo- 
nents v'g and v'^ are of particular interest because these 
may be compared with analogous correlations obtained 
from local helioseismology. Such correlations not only 
represent a potential diagnostic for giant-cell convec- 
tion, but they also reflect latitudinal angular momen- 



tum transport by Reynolds stresses which plays an es- 
sential role in maintaining the differential rotation pro- 
file ([21). However, the 2-D pdfs in Figure fTBb. d, and g 
appear nearly isotropic, implying that the horizontal ve- 
locity components near the surface {r — 0.98i?) are only 
weakly correlated. At high latitudes, there is a weak 
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Fig. 16. — Correlations between v'g-v'^ (left column), A' (middle column), and v'^—T' (right column) are shown at r = 0.98. (a)-(i) 
2-D pdfs averaged over a time interval of 62 days (yellow denotes the maximum contour level). Results are shown for latitude bands of 
0°-30° (g—i), 30°— 60° (d-f), and 60°— 90° (a— c), averaged over the northern and southern hemispheres (reversing the sign of v'g and in 
the southern hemisphere), (j)— (0 mean correlations, obtained by vertically averaging the 2-D histograms (pdfs) in each horizontal bin. 
Red, black, and blue lines represent low (0°-30°), mid (30°-60°), and high (60°-90°) latitudes respectively. 



positive correlation signifying equatorward angular mo- 
mentum transport, but at mid-latitudes the sense of the 
correlation reverses fFig. \W\i ). At low latitudes there is 
no clear systematic behavior, as expected if horizontal 



velocity correlations are induced by the vertical compo- 
nent of the rotation vector. 

The lack of prominent horizontal velocity correlations 
in the near-surface downflow network may be attributed 
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to the relatively small spatial and temporal scales of the 
convection. The effective Rossby number here is greater 
than for the larger-scale motions deeper in the convec- 
tion zone, implying weaker rotational influence (i i2.2p . 
Coriolis-induced correlations are consequently weaker. 
Since the differential rotation is maintained primarily 
by horizontal Reynolds stresses O, weaker horizontal 
velocity correlations help account for the decrease in lat- 
itudinal shear found in our simulation near the outer 
boundary. A near-surface shear layer is also found in 
helioseismic inversions, but in that case there is nearly 
uniform acceleration of the angular velocity at all lat- 
itudes so the latitudinal shea r does not change signifi- 
cantly (|Thompson et al.1l2003D . As discussed in 21 mean 
flows in the uppermost portion of the convection zone 
are likely sensitive to the subtle dynamics of the sur- 
face boundary layer and may require more sophisticated 
modeling approaches to capture fully. If a decrease in 
horizontal velocity correlations does indeed occur near 
the surface of the Sun as suggested by our simulations, 
then such correlations may be difficult to detect in helio- 
seismic inversions. 

A more promising diagnostic to search for in SSW 
maps may be correlations between horizontal divergence 
and vertical vorticity. Although there is much scat- 
ter, Figure [16] {b, e, h) demonstates a clear correlation 
between cyclonic vorticity and horizontal convergence 
which becomes more prominent at higher latitudes. This 
correlation is a signature of the Coriolis force as described 
in Sj3l There is also a correlation between anticyclonic 
vorticity and horizontal divergence but this only occurs 
for small values of The most intense vorticity of both 
signs occurs in regions of horizontal convergence. This 
is consistent with the interpretation discussed in S|3] in 
which downflow lanes are generally more turbulent than 
upflows. Vorticity of all orientations is generated by 
shear and entrainment and amplified by vortex stretch- 
ing^ 

'K omm et al.l ()2007t ) have presented evidence for cor- 
relations between C and A' in SSW maps derived from 
SOHO/MDI and GONG data. These correlations are 
approximately linear in regions of low magnetic activity, 
with cyclonic and anti-cyclonic vorticity associated with 
horizontal convergence and divergence respectively. This 
is qualitatively consistent with our simulation results 
since the high-amplitude anticyclonic vorticity in our 
simulations is associated with localized features which 
would be filtered out by the spatial averaging inherent in 
the helioseismic inversions. A more detailed comparison 
between our simulation results and SSW maps will be 
carried out in a subsequent paper. 

8. SUMMARY AND CONCLUSIONS 

High-resolution simulations of turbulent convection 
provide essential insight into the nature of global-scale 
motions in the solar convection zone, often referred to as 
giant cells, and into how these motions maintain the solar 
differential rotation and meridional circulation. Such in- 
sight is essential to inspire and interpret investigations of 
solar interior dynamics based on helioseismic inversions 
and photospheric observations. Although the simulation 
we focus on here is non-magnetic, our results have im- 
portant implications for solar dynamo theory and may 
be used to assess, calibrate, and further develop other 



modeling strategies such as mean-field models of solar 
and stellar activity cycles. 

The convective patterns realized in our simulations are 
intricate and continually evolving. Near the top of our 
computational domain at r = 0.98i?, there is an inter- 
connected network of downflow lanes reminiscent of pho- 
tospheric granulation but on a much larger scale. The 
power spectrum of the radial velocity peaks at ^ ~ 80, 
corresponding to a horizontal scale of about 50 Mm. 
However, a visual inspection of the convective patterns 
(Figs. [3 O in El) reveals a wide range of scales, with 
many cells spanning 10°-20° (100-200 Mm). Characteris- 
tic horizontal velocity scales are 250 m s~^ at r = 0.98i?, 
dropping to ^ 100 m in the mid convection zone. 
Near the surface, zonal flow amplitudes (v^) are on av- 
erage about 10% larger than latitudinal flow amplitudes 
(v'g) but in the mid convection zone all three velocity 
components have a comparable amplitude. Deep in the 
convection zone the surface network fragments into dis- 
connected downflow lanes and plumes but the skewness 
of the radial velocity remains strongly negative (Fig. 

A close inspection of the downflow network near the 
surface reveals a distinct tendency for structures to align 
in a north-south orientation at low latitudes. Such NS 
downflow lanes represent the largest and longest-lived 
features in the convection zone. Whereas correlation 
time scales for the downflow network are only a few days, 
NS downflow lanes can persist for weeks or even months. 
They are traveling convective modes which propagate in 
longitude about 8% faster than the equatorial rotation 
rate. Near the bottom of the shell the lanes fragment into 
downwelling plumes but some coherence extends across 
the entire convection zone (e.g. Fig.|4|). 

At higher latitudes, the downflow network is more 
isotropic and possesses intense cyclonic vorticity, induced 
by Coriolis forces. Localized cyclonic vortices are preva- 
lent near the interstices of the network at latitudes above 
about ±30°. These structures a re similar to the turbu- 
lent helical plumes observed by iBrummell et al.l ()1996f ) 
in Cartesian /-plane simulations and are associated with 
downward flow, horizontal convergence, and cool tem- 
peratures as well as cyclonic radial vorticity. They are 
confined to the upper convection zone (r > 0.95i?) and 
their horizontal scale is comparable to that of supergran- 
ulation, about 10-30 Mm. Typical lifetimes are several 
days to a week (e.g. Fig. [TT|) . High-latitude cyclones 
are highly intermittent and give rise to prominent expo- 
nential tails in the radial vorticity and temperature pdfs 

(Fig. [a. 

Near the surface, the horizontal divergence A' is highly 
correlated with the radial velocity v'^.. Thus, horizontal 
divergence fields obtained from SSW maps should pro- 
vide a good proxy for the radial velocity, at least on large 
scales. Furthermore, there is a strong correlation be- 
tween A' and the radial vorticity (' , with intense cyclonic 
vorticity in regions of horizontal convergence (downflows) 
amid a background of weaker anticyclonic vorticity in 
broader regions of divergence (upflows). This correla- 
tion applies over most of the horizontal surface area but 
breaks down for localized, high-amplitude events; the 
most intense vorticies, both cyclonic and anticyclonic, 
occur in downflow lanes (A' < 0). Correlations be- 
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tween C and A' represent a promising diagnostic for the 
investigat ion of lar ge-scale flow patterns in SSW maps 
( |Komm ct~aLll2007[) . 

A significant new feature of the simulation presented 
here relative to previous models is the manner in which 
the differential rotation is maintained. As demonstrated 
in Figure [71 the resolved convective motions transport 
angular momentum equatorward and inward by means 
of Reynolds stresses while the meridional circulation op- 
poses this transport, such that 

v.(FR^TFMcy « . (18) 

This simulation has thus crossed a threshold in which 
viscous diffusion no longer contributes significantly 
to the angular momentum balance. 

The implications of this result are profound. Although 
there are subtle nonlinear feedbacks, the meridional cir- 
culation pattern is largely determined by F^'^ under the 
constraint that the resulting mean flows satisfy equa- 
tion Convective motions (particularly NS down- 
flow lanes) redistribute angular momentum and the re- 
sulting differential rotation induces circulations through 
Coriolis forces until a steady state is reached. Baroclin- 
icity also plays an important role, breaking the Taylor- 
Proudman constraint which favors cylindrical rotation 
profiles. Baroclinic torques arise in part from thermal 
coupling to the tachocline which is represented in our 
simulation by imposing a latitudinal entropy gradient on 
the lower boundary at the base of the convection zone 
(Q. 

The mean flows which result are similar to those in- 
ferred from helioseismic inversions. The mean angu- 
lar velocity decreases monotonically with latitude with 
nearly radial contours at mid-latitudes. This is similar 
to the solar rotation profile, although the angular ve- 
locity contrast of ~ 50 nHz be tween 0°-60° is smalle r 
than the - 90 nHz in the Sun (jThompson et al.ll2003h . 
The time-averaged meridional circulation is dominated 
by a single cell in each hemisphere with poleward flow of 
about 20 m in the upper convection zone (notwith- 
standing a flow reversal near the upper boundary which 
may be artificial). The sense and amplitude of this cir- 



culation is comparable to that inferred near the sur- 
face of the Sun from Doppler measurements and helio- 
seismic inversions. The lower convection zone currently 
lies beyond the reach of helioseismic probing but many 
have proposed that an equatorward return fiow may ex- 
ist and furthermore that this global circulation pattern 
may play an esse ntial role in est ablishing th e solar ac- 
tivity cycle (Dikpati fc Charbonneau 1999; Di kpati et al.l 
12004 e.g.). For a revi ew of these so- c alled flux-transport 
dynamo models, see ICharbonneaul (|2005f l. The mean 
meridional circulation in our simulation is similar to that 
used in many flux-transport dynamo models but month- 
to-month fluctuations about this mean are large. 

In order to gain further insight into how the global 
solar dynamo operates and into how mean flows are 
maintained, we must extend the lower boundary of our 
computational domain below the solar convection zone 
and thus explicitly resolve the complex dynamics oc- 
curring in the overshoot region and tachocline. Fur- 
thermore, we must incorporate magnetism and investi- 
gate how magnetic flux is amplified, advected, and or- 
ganized by turbulent penetrative convection, rotational 
shear, and global circulations. S uch efforts are already 
underway (jBrowning et al.ll2006l ) and will continue. Fu- 
ture work will also focus on improving our understand- 
ing of the upper convection zone, including how granula- 
tion and supergranulation influence giant cells and mean 
flow patterns, and how signatures of internal flows and 
magnetism might be manifested in helioseismic measure- 
ments and photospheric observations. 
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